
# ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
# DESCRIPTION
# ______________________________________________________________________________

# Author: Changwook Ju, Yale University (changwook.ju@yale.edu)

# This code creates Figures 2 and 3.

# Preamble ---------------------------------------------------------------------

rm(list = ls())

# Packages ---------------------------------------------------------------------

library(tidyverse)
library(ggpmisc)
library(ggpubr)
library(grid)
library(haven)
library(patchwork)
library(signs)

# Plot -------------------------------------------------------------------------

# edit geomsegment legend key

GeomSegment$draw_key <-  function (data, params, linewidth) {
  draw_key_vpath <- function (data, params, linewidth) {
    segmentsGrob(0.1,
                 0.5,
                 0.9,
                 0.5,
                 gp = gpar(
                   col = alpha(data$colour, data$alpha),
                   lwd = data$linewidth * .pt,
                   lty = data$linetype
                 ))
  }
  grobTree(draw_key_vpath(data, params, linewidth),
           draw_key_point(data, params, linewidth))
}

## Child Soldiering: US DoS ----------------------------------------------------

# data

cs_usdos_ol <-
  read_dta("./results/cs_usdos_ol.dta") %>%
  filter(!parm %in% c("cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(eq == "statedept_sv" ~ "OL"))

cs_usdos_ziol <-
  read_dta("./results/cs_usdos_ziol.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    model_name = case_when(
      eq == "statedept_sv" ~ "ZIOL (Outcome)",
      eq == "inflate" ~ "ZIOL (Inflation)"
    )
  )

cs_usdos <-
  rbind(cs_usdos_ol, cs_usdos_ziol)

cs_usdos <- cs_usdos %>%
  mutate(
    new_parm = case_when(
      parm == "csindex" ~ "Child Soldiering",
      parm == "secession" ~ "Secessionist",
      parm == "islam" ~ "Islamist",
      parm == "com" ~ "Communist",
      parm == "central_control" ~ "Central Control",
      parm == "territory" ~ "Territorial Control",
      parm == "state_support" ~ "Foreign Support",
      parm == "lootable1" ~ "Lootable Resources",
      parm == "duration" ~ "Duration",
      parm == "autocracy" ~ "Regime Type",
      parm == "side_a_sexviolence" ~ "Government\nSexual Violence"
    ),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

cs_usdos$new_parm <-
  factor(cs_usdos$new_parm,
         levels = rev(unique(cs_usdos$new_parm)))

cs_usdos$model_name <-
  factor(cs_usdos$model_name,
         levels = rev(unique(cs_usdos$model_name)))

# plot

space_between_bars <- 0.15
cs_usdos$type_y_adj =
  scale(as.numeric(as.factor(cs_usdos$model_name))) * space_between_bars
cs_usdos$y =
  as.numeric(as.numeric(as.factor(cs_usdos$new_parm)) + cs_usdos$type_y_adj)

cs_usdos_plot <- cs_usdos %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (US DoS)\nObservations = 299") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:11,
    labels = rev(unique(cs_usdos$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) + coord_cartesian(clip = "off")

cs_usdos_plot_inset <- cs_usdos %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  scale_x_continuous(
    breaks = seq(-0.3, 0.3, 0.3),
    limits = symmetric_limits,
    labels = signs_format()
  ) +
  scale_y_continuous(
    breaks = 1:11,
    labels = rev(unique(cs_usdos$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("p < 0.05"),
               "p \u2265 0.05"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    plot.background = element_blank()
  ) +
  coord_cartesian(xlim = c(-0.3, 0.3),
                  ylim = c(3, 3))

cs_usdos_plot <- cs_usdos_plot +
  annotate(
    "segment",
    x = 4.35,
    xend = 0.75,
    y = 3,
    yend = 3,
    color = "black",
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.75,
    xmax = 0.75,
    ymin = 2.5,
    ymax = 3.5,
    fill = NA,
    color = "black",
    linewidth = 0.2
  ) +
  annotation_custom(
    grob = ggplotGrob(cs_usdos_plot_inset),
    xmin = 2,
    xmax = 14,
    ymin = 1.95,
    ymax = 3.8
  )

## Forced Child Recruitment: US DoS --------------------------------------------

# data

fcr_usdos_ol <-
  read_dta("./results/fcr_usdos_ol.dta") %>%
  filter(!parm %in% c("cut1", "cut2", "cut3")) %>%
  mutate(model_name = case_when(eq == "statedept_sv" ~ "OL"))

fcr_usdos_ziol <-
  read_dta("./results/fcr_usdos_ziol.dta") %>%
  filter(!parm %in% c("_cons", "cut1", "cut2", "cut3")) %>%
  mutate(
    model_name = case_when(
      eq == "statedept_sv" ~ "ZIOL (Outcome)",
      eq == "inflate" ~ "ZIOL (Inflation)"
    )
  )

fcr_usdos <-
  rbind(fcr_usdos_ol, fcr_usdos_ziol)

fcr_usdos <- fcr_usdos %>%
  mutate(
    new_parm = case_when(
      parm == "forceindex" ~ "Forced Child Recruitment",
      parm == "secession" ~ "Secessionist",
      parm == "islam" ~ "Islamist",
      parm == "com" ~ "Communist",
      parm == "central_control" ~ "Central Control",
      parm == "territory" ~ "Territorial Control",
      parm == "state_support" ~ "Foreign Support",
      parm == "lootable1" ~ "Lootable Resources",
      parm == "duration" ~ "Duration",
      parm == "autocracy" ~ "Regime Type",
      parm == "side_a_sexviolence" ~ "Government\nSexual Violence"
    ),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

fcr_usdos$new_parm <-
  factor(fcr_usdos$new_parm,
         levels = rev(unique(fcr_usdos$new_parm)))

fcr_usdos$model_name <-
  factor(fcr_usdos$model_name,
         levels = rev(unique(fcr_usdos$model_name)))

# plot

space_between_bars <- 0.15
fcr_usdos$type_y_adj =
  scale(as.numeric(as.factor(fcr_usdos$model_name))) * space_between_bars
fcr_usdos$y =
  as.numeric(as.numeric(as.factor(fcr_usdos$new_parm)) + fcr_usdos$type_y_adj)

fcr_usdos_plot <- fcr_usdos %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (US DoS)\nObservations = 226") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:11,
    labels = rev(unique(fcr_usdos$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

fcr_usdos_plot_inset <- fcr_usdos %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 0.75,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 2.5,
    stroke = 1,
    fill = "white"
  ) +
  scale_x_continuous(
    breaks = seq(-0.4, 0.4, 0.4),
    limits = symmetric_limits,
    labels = signs_format()
  ) +
  scale_y_continuous(
    breaks = 1:11,
    labels = rev(unique(fcr_usdos$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("p < 0.05"),
               "p \u2265 0.05"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    axis.title = element_blank(),
    axis.ticks.y = element_blank(),
    axis.ticks.x = element_line(color = "black"),
    axis.text.y = element_blank(),
    axis.text.x = element_text(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    plot.background = element_blank()
  ) +
  coord_cartesian(xlim = c(-0.4, 0.4),
                  ylim = c(3, 3))

fcr_usdos_plot <- fcr_usdos_plot +
  annotate(
    "segment",
    x = 2.375,
    xend = 0.5,
    y = 3,
    yend = 3,
    color = "black",
    linewidth = 0.2
  ) +
  annotate(
    "rect",
    xmin = -0.5,
    xmax = 0.5,
    ymin = 2.5,
    ymax = 3.5,
    fill = NA,
    color = "black",
    linewidth = 0.2
  ) +
  annotation_custom(
    grob = ggplotGrob(fcr_usdos_plot_inset),
    xmin = 1.3,
    xmax = 9.15,
    ymin = 1.95,
    ymax = 3.8
  )

## Child Soldiering: AI --------------------------------------------------------

# data

cs_ai_ol <-
  read_dta("./results/cs_ai_ol.dta") %>%
  filter(parm %in% "csindex") %>%
  mutate(model_name = case_when(eq == "ai_sv" ~ "OL"))

cs_ai_ziol <-
  read_dta("./results/cs_ai_ziol.dta") %>%
  filter(parm %in% "csindex") %>%
  mutate(model_name = case_when(
    eq == "ai_sv" ~ "ZIOL (Outcome)",
    eq == "inflate" ~ "ZIOL (Inflation)"
  ))

cs_ai <-
  rbind(cs_ai_ol, cs_ai_ziol)

cs_ai <- cs_ai %>%
  mutate(
    new_parm = case_when(parm == "csindex" ~ "Child Soldiering"),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

cs_ai$new_parm <-
  factor(cs_ai$new_parm,
         levels = rev(unique(cs_ai$new_parm)))

cs_ai$model_name <-
  factor(cs_ai$model_name,
         levels = rev(unique(cs_ai$model_name)))

# plot

space_between_bars <- 0.15
cs_ai$type_y_adj =
  scale(as.numeric(as.factor(cs_ai$model_name))) * space_between_bars
cs_ai$y =
  as.numeric(as.numeric(as.factor(cs_ai$new_parm)) + cs_ai$type_y_adj)

cs_ai_plot <- cs_ai %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (AI)\nObservations = 299") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:1,
    labels = rev(unique(cs_ai$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

cs_ai_plot

## Forced Child Recruitment: AI ------------------------------------------------

# data

fcr_ai_ol <-
  read_dta("./results/fcr_ai_ol.dta") %>%
  filter(parm %in% "forceindex") %>%
  mutate(model_name = case_when(eq == "ai_sv" ~ "OL"))

fcr_ai_ziol <-
  read_dta("./results/fcr_ai_ziol.dta") %>%
  filter(parm %in% "forceindex") %>%
  mutate(model_name = case_when(
    eq == "ai_sv" ~ "ZIOL (Outcome)",
    eq == "inflate" ~ "ZIOL (Inflation)"
  ))

fcr_ai <-
  rbind(fcr_ai_ol, fcr_ai_ziol)

fcr_ai <- fcr_ai %>%
  mutate(
    new_parm = case_when(parm == "forceindex" ~ "Forced Child Recruitment"),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

fcr_ai$new_parm <-
  factor(fcr_ai$new_parm,
         levels = rev(unique(fcr_ai$new_parm)))

fcr_ai$model_name <-
  factor(fcr_ai$model_name,
         levels = rev(unique(fcr_ai$model_name)))

# plot

space_between_bars <- 0.15
fcr_ai$type_y_adj =
  scale(as.numeric(as.factor(fcr_ai$model_name))) * space_between_bars
fcr_ai$y =
  as.numeric(as.numeric(as.factor(fcr_ai$new_parm)) + fcr_ai$type_y_adj)

fcr_ai_plot <- fcr_ai %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (AI)\nObservations = 226") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:1,
    labels = rev(unique(fcr_ai$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

fcr_ai_plot

## Child Soldiering: Max -------------------------------------------------------

# data

cs_max_ol <-
  read_dta("./results/cs_max_ol.dta") %>%
  filter(parm %in% "csindex") %>%
  mutate(model_name = case_when(eq == "max_sv" ~ "OL"))

cs_max_ziol <-
  read_dta("./results/cs_max_ziol.dta") %>%
  filter(parm %in% "csindex") %>%
  mutate(model_name = case_when(
    eq == "max_sv" ~ "ZIOL (Outcome)",
    eq == "inflate" ~ "ZIOL (Inflation)"
  ))

cs_max <-
  rbind(cs_max_ol, cs_max_ziol)

cs_max <- cs_max %>%
  mutate(
    new_parm = case_when(parm == "csindex" ~ "Child Soldiering"),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

cs_max$new_parm <-
  factor(cs_max$new_parm,
         levels = rev(unique(cs_max$new_parm)))

cs_max$model_name <-
  factor(cs_max$model_name,
         levels = rev(unique(cs_max$model_name)))

# plot

space_between_bars <- 0.15
cs_max$type_y_adj =
  scale(as.numeric(as.factor(cs_max$model_name))) * space_between_bars
cs_max$y =
  as.numeric(as.numeric(as.factor(cs_max$new_parm)) + cs_max$type_y_adj)

cs_max_plot <- cs_max %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (Max)\nObservations = 299") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:1,
    labels = rev(unique(cs_max$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

cs_max_plot

## Forced Child Recruitment: Max -----------------------------------------------

# data

fcr_max_ol <-
  read_dta("./results/fcr_max_ol.dta") %>%
  filter(parm %in% "forceindex") %>%
  mutate(model_name = case_when(eq == "max_sv" ~ "OL"))

fcr_max_ziol <-
  read_dta("./results/fcr_max_ziol.dta") %>%
  filter(parm %in% "forceindex") %>%
  mutate(model_name = case_when(
    eq == "max_sv" ~ "ZIOL (Outcome)",
    eq == "inflate" ~ "ZIOL (Inflation)"
  ))

fcr_max <-
  rbind(fcr_max_ol, fcr_max_ziol)

fcr_max <- fcr_max %>%
  mutate(
    new_parm = case_when(parm == "forceindex" ~ "Forced Child Recruitment"),
    sig = ifelse(p < 0.05, 1, 0.2)
  )

fcr_max$new_parm <-
  factor(fcr_max$new_parm,
         levels = rev(unique(fcr_max$new_parm)))

fcr_max$model_name <-
  factor(fcr_max$model_name,
         levels = rev(unique(fcr_max$model_name)))

# plot

space_between_bars <- 0.15
fcr_max$type_y_adj =
  scale(as.numeric(as.factor(fcr_max$model_name))) * space_between_bars
fcr_max$y =
  as.numeric(as.numeric(as.factor(fcr_max$new_parm)) + fcr_max$type_y_adj)

fcr_max_plot <- fcr_max %>%
  ggplot(aes(y = new_parm)) +
  geom_vline(xintercept = 0,
             linetype = 2) +
  geom_segment(
    aes(
      x = min95,
      xend = max95,
      y = y,
      yend = y,
      color = model_name,
      alpha = sig
    ),
    linewidth = 1,
    lineend = "round"
  ) +
  geom_point(
    aes(
      x = estimate,
      y = y,
      color = model_name,
      alpha = sig,
      shape = model_name
    ),
    size = 3,
    stroke = 1.5,
    fill = "white"
  ) +
  guides(
    color = guide_legend(reverse = TRUE, order = 1),
    shape = guide_legend(reverse = TRUE, order = 1),
    alpha = guide_legend(override.aes = list(shape = c(NA, NA)), order = 2)
  ) +
  facet_wrap(paste0("Rebel CRSV (Max)\nObservations = 226") ~ .,
             scales = "free_y") +
  labs(color = "Model (Stage)",
       shape = "Model (Stage)",
       alpha = paste("Significance")) +
  scale_x_continuous(limits = symmetric_limits,
                     labels = signs_format()) +
  scale_y_continuous(
    breaks = 1:1,
    labels = rev(unique(fcr_max$new_parm)),
    expand = expansion(add = 0.5)
  ) +
  scale_alpha_identity(
    breaks = c(1, 0.2),
    labels = c(paste("Yes (p < 0.05)"),
               "No (p \u2265 0.05)"),
    guide = "legend"
  ) +
  scale_shape_manual(values = c(23,
                                22,
                                21)) +
  scale_color_manual(values = c("blue",
                                "red",
                                "limegreen")) +
  theme_bw() +
  theme(
    strip.placement = "outside",
    strip.background = element_rect(color = "black", fill = "white"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black"),
    panel.border = element_rect(fill = NA, color = "black"),
    panel.grid.minor = element_blank(),
    axis.title = element_blank(),
    axis.ticks.y  = element_blank(),
    legend.position = "none",
    plot.margin = margin(10, 10, 10, 10, "pt"),
    legend.background = element_blank(),
    legend.box.background = element_blank(),
    legend.margin = margin(0, 3, 0, 3, "pt"),
    legend.spacing.x = unit(3, "pt"),
    legend.key.width = unit(30, "pt"),
    legend.box.margin = margin(0, 0, 0, 0, "pt"),
    strip.text.x = element_text(margin = margin(7, 0, 7, 0, "pt"),
                                lineheight = 1.2)
  ) +
  coord_cartesian(clip = "off")

fcr_max_plot

# Patchwork --------------------------------------------------------------------

combined_1 <-
  (cs_usdos_plot | fcr_usdos_plot) &
  theme(
    legend.position = "bottom",
    plot.margin = margin(5, 5, 5, 5, "pt"),
    text = element_text(size = 15)
  )

combined_1 + plot_layout(guides = "collect")

ggsave(
  "figures/coefficient plots/coefficient_plot_1.pdf",
  device = cairo_pdf,
  width = 12,
  height = 10
)

combined_2 <-
  ((cs_ai_plot | fcr_ai_plot) / (cs_max_plot | fcr_max_plot)) &
  theme(
    legend.position = "bottom",
    plot.margin = margin(5, 5, 5, 5, "pt"),
    text = element_text(size = 15)
  )

combined_2 + plot_layout(guides = "collect")

ggsave(
  "figures/coefficient plots/coefficient_plot_2.pdf",
  device = cairo_pdf,
  width = 12,
  height = 4.75
)
